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Abstract. A second-order Taylor series expansion of the free energy functional 
provides analytical expressions for the magnetic field dependence of the free energy 
and of the magnetization of ferrofluids, here modelled by dipolar Yukawa interaction 
potentials. The corresponding hard core dipolar Yukawa reference fluid is studied 
within the framework of the mean spherical approximation. Our findings for 
the magnetic and phase equilibrium properties are in quantitative agreement with 
previously published and new Monte Carlo simulation data. 



1. Introduction 

Ferrofluids are colloidal suspensions of ferromagnetic grains dispersed in a solvent. 
To keep such dispersions stable, one must avoid aggregation and counterbalance the 
attractive van der Waals interaction by repulsive forces, either electrostatic ones, if 
the particles are charged, or steric ones, if their surface is coated with polymers or 
surfactants. Each particle of a ferrofluid possesses a permanent magnetic dipole moment 
so that the ferrocolloids have much in common with dipolar liquids [HIS]- The effective 
interaction of such magnetic particles can be modeled by the Stockmayer potential, 
which is the sum of an isotropic Lennard- Jones (LJ) and the anisotropic dipole-dipole 
interaction [3l H] . Within this approach the L J interaction parameters are considered to 
be effective ones, incorporating the effects of the solvent. (In a more reflned approach 
one would model the suspension as a binary mixture of a nonpolar solvent and a polar 
solute component [5l[6].) 

Based on previous experience [3, E] , here we substitute the Lennard- Jones potential 
by a hard core Yukawa potential. For a homogeneous isotropic dipolar Yukawa (DY) 
fluid this substitution provides the valuable beneflt that one can solve analytically 
the corresponding Ornstein-Zernike equation within the mean spherical approximation 
(MSA). It turns out that the thermodynamic properties of a DY fluid can be calculated 
easily from the corresponding results for a hard core Yukawa and a dipolar hard sphere 



Magnetization and susceptibility of ferrofluids 



2 



(DHS) fluid. The thermodynamic properties of a DY fluid differ from those of the DHS 
fluid but the magnetic susceptibility of the DY fluid is the same as that of the DHS 
fluid, which was predicted first by Wertheim [9j. The disadvantage of the MSA is that 
it can predict the magnetization only for weak magnetic fields because it is a linear 
response theory. In order to eliminate this shortcoming, first Lebedev [10] suggested 
a semi-empirical method for using MSA at finite magnetic fields. Later Morozov and 
Lebedev [IT] extended the applicability of DHS MSA to the case of arbitrary strengths of 
the magnetic fields on the basis of the Lovett-Mou-Buff equation [12], relating the single 
particle distribution function and the direct correlation function. Within the framework 
of density functional theory, here we also propose an equation for the magnetization of 
ferrofluids, which turns out to be simpler than the corresponding equation of Morozov 
and Lebedev [TT] . In order to facilitate the calculation of the phase diagrams, we present 
also the magnetic field dependence of the chemical potential and of the pressure. 



2. The model 



We study so-called hard core dipolar Yukawa fluids which consist of spherically 
symmetric particles interacting via a Yukawa (Y) potential characterized by parameters 
cr, e, and A: 

uriru) = -— exp[-A(ri2 - a)]. (1) 

f 12 

In addition there is a dipolar interaction due to point dipoles embedded at the centres 
of the particles: 

2 

Tfl 

MD(ri2, UJi,UJ2) = —D{uJi2, LOi, UJ2)), (2) 

with the rotationally invariant function 

D{uJi2,uJi,uj2) = 3[mi(u;i) ■ ris] [in2(cu2) -^12] - [mi(a;i) -£12(6^2)], (3) 

where particle 1 (2) is located at ri (r2) and carries a dipole moment of strength m with 
an orientation given by the unit vector mi(co'i) (m2(co'2)) with polar angles Ui = (6*1, 0i) 
(0^2 = (6'2,02)); i'i2 = i"i — r2 is the difference vector between the centres of particle 1 
and 2, ri2 = |ri2|, and ru = Tu/r^ is a unit vector with orientation luu = (6'i2,0i2)- 
The hard core dipolar Yukawa potential is defined by the aforementioned interaction 
potentials as 

J 00 , r < (J 

UDY{ri2,UJi,UJ2) = < / N , / N ^ 4 

[ UY{ri2) +UD{ri2,uJi,uJ2) , r > cr. 

In the following we consider the DY fluid in a homogeneous external magnetic field 
H (the direction of which is taken to coincide with the direction of the z axis). The 
contribution of the latter to the interaction potential is 

Uexti^) = —mHcosO, (5) 

where the angle 6 measures the orientation of the dipole of a single particle relative to 
the field direction. 



Magnetization and susceptibility of ferrofluids 



3 



3. Homogeneous magnetization in an external field 

In general, an inhomogeneous and anisotropic dipolar fluid is described by the one 
particle distribution function p{r,uj). For bulk fluid phases the number density p = 
J dojpijc, oj) is spatially constant and Q;(r, uj) = p(r, u;)/p is the orientational distribution 
function at the point r. In a homogeneously magnetized bulk phase the orientational 
distribution function depends only on u {a{r,uj) = C({uj)) and it is normalized to 1, i.e., 
J duja{uj) = 1. In these terms our subsequent bulk analysis is based on the minimization 
of the following grand canonical variational functional: 

^l[p,{a{uj)},T, p] = FDY[p,Wi^)}^T] - J d^rdujpa{uj){p- Uextiuj)), (6) 

where Fdy is the Helmholtz free energy functional of the DY fluid, T is the temperature 
and p is the chemical potential. The Helmholtz free energy functional consists of an 
ideal gas and an excess contribution: 

Fnrip, {«H},T] = F"^[p, {aM},T] + F^Jy[p, {a{u;)},T]. (7) 

For the system under study the ideal contribution has the form 

F''^[p, {a{uj)},T] = Vp-^p ln(pA3) + J dua{u) ln(47ra(cu)) , (8) 

where A is the de Broglie wavelength, /3 = l/{kBT) is the inverse temperature and 
V is the volume of the fluid. The central quantity in this theory is the excess free 
energy functional F*^^, which originates from interparticle interactions. The excess DY 
free energy functional for an anisotropic system is not known. However, it can be 
approximated by a functional Taylor series, expanded around a homogeneous zsotropic 
reference system with bulk density p. Neglecting all terms beyond second order one has 

/?F^r[P,{«H}>2^] = /?^Dy(P,T)-p I d'r,dcoAaico)c^^Up) 

J d^ridui J d^r2du2Aa{ui)Aa{u2)c^^Y{ri2,uJi,U2, p), (9) 

where Aa^u) = a{uj) — l/(47r) is the difference between the anisotropic and the 
isotropic orientational distribution function. The flrst- (c^^^) and second-order (c*-^-') 
direct correlation functions appearing in Eq. iQ are the flrst and second functional 
derivatives of [—(5F^y) evaluated at the homogeneous and isotropic density p. Since 
j dujAa{uj) = only the second- and higher-order direct correlation functions provide 
nonzero contributions to the functional j3F'j^Y- the following we approximate the 
second-order direct correlation function of a DY fluid by the corresponding MSA solution 
due to Henderson a/ [7]: 

('2') ('2') ('2) ('2) 

C})y(ri2,C^l,C^2,p) = Cy{ri2,p) + Cy {ri2, p)D{uJi2,UJi,UJ2) + c\'{ri2, p)A{uJi,UJ2) (10) 

(see Eq. (jS])) with the rotationally invariant function 

A{uJi,uj2) = rhi{uji) ■m2{uj2)- (11) 



p! 
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In Eq. ( fTOj) . where Cy is the hard core Yukawa direct correlation function [13], c}^ and 

(2) 

c]^ are spherical harmonic components of the corresponding direct correlation function 
of a DHS fluid in MSA [9]. The latter two functions can be expressed in terms of the 
Percus-Yevic hard sphere direct correlation function. If the volume of the fluid V has 
the shape of a rotational ellipsoid (elongated around the magnetic field direction) and 
the magnetization is homogeneous, a{uj) depends only on the angle 6 relative to the long 
axis of the ellipsoid, and thus it can be expanded in terms of Legendre polynomials: 



a{uj) = ^a{cose) = ^J^aiP^{cose). (12) 

Due to ao = 1/2 one has Aa{uj) = ^ ctiPi^cos 9) . In order to avoid depolarization 
effects, we consider sample shapes of infinitely prolate eUipsoids, i.e., needle-shaped 
volumes V. Since in an expansion in terms of Pi{cosd) the orthogonal functions D and 
A (see Eq. (ITOl) ) have contributions only of the order i < 1, ai is the only term in 
a(cos6') which contributes to Eq. Thus after some elementary calculations one 

finds for the free energy functional of Eq. iQ 

fS>.W}.r| ,^„(^_r)-|(l-,(-0)a;. (13) 

where /|,y = F^y/^ is the Helmholtz free energy density of the isotropic DY fluid and 
. N (l + 2x)2 ^^^^ 
(1 — x)^ 

is the reduced inverse compressibility function of hard spheres within the Percus-Yevic 
approximation [H]. The parameter ^(xl) stems from the DHS MSA and is given by 
the implicit equation 

A7rxL = q{20-q{-0, (15) 

where xl = {Spiv? /?> is the Langevin susceptibility. Minimization of the grand canonical 
functional with respect to the orientational distribution function yields 

a(cos0) = C-^ exp [{l3mH + 2(1 - g(-0)ai)Pi(cos^)] , (16) 

with the normalization constant 

^ 2sinh(/3mg + 2(l-g(-0)ai) 

j3mH + 2{l-q{-i))ai ' ^ ^ 

The coefficient ai is determined by the implicit equation 

ai = ^-L {(3mH + 2(1 - g(-0)«i) , (18) 

where L{x) is the well known Langevin function. Further details of this minimization 
scheme can be found in Ref. [1]. Minimization of Q/V with respect to p yields 

(3p = f3pDY{p, T) - ln[C/2] + + ^(1 - (19) 
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where fi dy is the chemical potential of the DY fluid [H [S] . By eliminating the chemical 
potential the equilibrium grand canonical potential can be expressed as 

y = /wlP, T) - pfiDYiP, T) - -—a^ , (20) 

where /dy = F^y /V is the free energy density of the homogeneous and isotropic DY 
fluid [71 [8] . Due to = —pV the corresponding pressure of an equilibrium phase as a 
function of the density and temperature is 

P = PDY{p,T) + -ja^ , (21) 

where Pdy = —Idy+PPdy is the pressure of the DY fluid. Each particle of the magnetic 
fluid carries a dipole moment which may be preferentially aligned in a particular 
direction. This gives rise to a magnetization M (the direction of which coincides with 
the direction of the external magnetic field): 

f 2 

M = p duja{uj)m cos 6 = -mpai. (22) 

Equations (l22i) and (ITSj) lead to an implicit equation for the external field dependence 
of the magnetization: 

M = mpL (pmH + SM^^^^^tl^ ) . (23) 



\ mp 
Accordingly, the zero field magnetic susceptibility x is 

This result is identical with the ones of Henderson et al [7] for DY fluids and with 
Wertheim's [9] equation for DHS fluids. 

4. Monte Ccirlo simulations 

We carried out Monte Carlo simulations for DY fluids using NVT and NpT ensembles 
in order to verify the predictions of the present DFT. The initial configurations in 
a cubic box were created by random insertions of the particles, excluding particle 
overlaps. Boltzmann sampling and periodic boundary conditions with the minimum- 
image convention were applied. A spherical cutoff of the pair interaction potentials 
at half of the cell length was applied, and long-ranged corrections (LRC) were taken 
into account. For the dispersion part of the potential, the usual Yukawa-tail LRC 
was used while for the dipole-dipole interaction a reaction field LRC and a conducting 
boundary condition were applied. For the simulation of the magnetization after 40000 
equilibration cycles 0.6-0.8 million production cycles were used involving 512 particles. 
In the simulations the equilibrium magnetization can be obtained from the expression 

M = ^/E-^)' (25) 
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where the brackets denotes the ensemble average. 

We obtained the external field dependence of the numerical data for liquid-vapor 
coexistence of a DY fluid by using the extended NpT plus test particle (NpT+TP) 
method, which is described in detail in Ref. [16] so that here we present only an outline. 
Starting from a point (To,po) in the one-phase region of the {T,p) parameter plane, the 
chemical potential can be expanded into a Taylor series about the point (To,po) up 
to third order: 

% • 

i=l 

On the basis of basic thermodynamic relations the coefficients of this series can be 
expressed in terms of the derivatives of the enthalpy and the volume of the system with 
respect to T and p and they can be calculated from fiuctuation formulas by performing 
an NpT+TP Monte Carlo simulation at the point {Tq,pq). These derivatives and 
fiuctuation formulas are given in Ref. [16]. Carrying out this scheme for a point in the 
liquid and in the vapor phase, respectively, and rewriting the third-order Taylor series of 
fi for these points accordingly, the vapor pressure curve as well as other equilibrium data 
can be obtained from the intersection of these curves in the appropriate temperature 
range within an accuracy depending on the number of terms taken into account in 
the Taylor expansion. The NpT ensemble MC simulations involved 512 particles and 
about 0.8 million cycles were performed. The chemical potential was calculated by using 
Widom's test particle method. 

5. Results and discussion 

In the following we shall use reduced quantities: T* = kT/e as reduced temperature, 
p* = pa^ as reduced density, m* = m/y/eo^ as reduced dipole moment, H* = H^/o^Je 
as reduced magnetic field strength, and M* = M^fo^t as the reduced magnetization. 
For the reference system, all results are obtained for A = 1.8/cr (see Eq. ([T])), and 
different from Ref. [8] the so-called first-order MSA [15] is used to describe the Yukawa 
fiuid. The magnetization curves obtained from the numerical solution of Eqs. fl25]) and 
( IT5|) for a needle-shaped sample are displayed in Fig. 1(a) together with the result of 
the Langevin theory for m* = 1. Figure 1(a) shows that the interparticle interaction 
enhances the magnetization relative to the corresponding values of the Langevin theory. 
The theoretical data of the present theory are consistent with our new simulation data. 
For the reduced density p* = 0.25 we find both DFT and MC results to be in excellent, 
quantitative agreement for all field strengths. For p* = 0.6 the level of quantitative 
agreement reduces for larger field strengths. 

The density dependence of the magnetic susceptibility obtained from the numerical 
solution of Eqs. fl23l) and f|T5l) in comparison with earlier simulation data of Szalai et 
al [8j is displayed in Fig. 1(b). For the reduced dipole moment m* = 1 the agreement 
between the theoretical and simulation results is quite good, especially at lower densities. 
For m* = the level of quantitative agreement reduces for larger reduced densities. 



(^-^o)^ + (p-Po)^ 



/i(T,p). (26) 
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Figure 1. (a) Magnetization curves of DY fluids for different values of the reduced 

density p* and reduced temperature T* for m* = 1. (b) Magnetic zero- field 
susceptibility of DY fluids for different values of the reduced dipole moment and for 
the reduced temperature T* = 1. The symbols in both flgures are the results of Monte 
Carlo simulations whereas the lines are obtained from the present DFT. 




Figure 2. (a) Evolution of the phase diagram as function of the external magnetic 
field for a DY fluid in a needle-shaped container for m* = 1. (b) Shift of the critical 
temperature (AT* = T*{H) — T*(0)) of the liquid- vapor critical point as a function 
of the applied field for m* = 1. The symbols in both figures are the results of Monte 
Carlo simulations whereas the lines are obtained from the present DFT. 



The liquid-vapor phase diagrams at a given temperature T follow from requiring 
the equality of the chemical potential and of the pressure for the densities of coexisting 
phases. The phase diagrams in the density-temperature plane are presented in Fig. 
2(a) for various fixed values of the magnetic field and for the reduced dipole moment 
m* — 1. Figure 2(a) shows that for H* — and m* — 1 there is no spontaneous 
magnetization within the framework of this theory, i.e., there is no ferromagnetic liquid 
phase. For H* there is a first-order phase transition between a low density gas 
and a high density liquid phase and both phases exhibit a nonzero magnetization. The 
agreement between the NpT-t-TP simulation and the DFT data is reasonable. For 
m* — 1 Fig. 2(b) represents the magnetic field dependence of the critical temperature 
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corresponding to the liquid-gas phase transitions, obtained from the present DFT and 
hquid-gas equihbrium NpT+TP simulation data. For small H* both the theoretical and 
the simulation shifts increase quadratically as a function of the field strength. Similar 
results have been found for Stockmayer fluids in Refs. [TTJ [18] . 

6. Summary 

The following main results have been obtained. 

(1) Based on a second-order Taylor series expansion of the free energy functional 
of the dipolar Yukawa fluids, within the mean spherical approximation we have derived 
an analytical expression for the magnetization of ferrofluids. There is good agreement 
between the results of the density functional theory and the Monte Carlo simulation 
data for reduced dipole moments m* < 1.4. The proposed implicit equation for the 
magnetization (Eq. ( l23l) ) extends the applicabihty of MSA for arbitrary strengths of 
the magnetic fields. 

(2) We have derived analytical equations for the external field dependence of the 
pressure (Eq. fl2T|) ) and of the chemical potential (Eq. f|T9l) ) for dipolar Yukawa fluids 
which allows one to determine the magnetic field dependence of the liquid-vapor phase 
equilibria. 
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